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Abstract 

The analytical gradient for periodic systems is presented, for the case of 
metallic systems. The total energy and the free energy are computed on the 
Hartree-Fock or density functional level, with the wave function being ex- 
panded in terms of Gaussian type orbitals. The expression for the gradient 
is similar to the case of insulating systems, when no thermal broadening is ap- 
plied. When the occupation of the states is according to the Fermi function, 
then the gradient is consistent with the gradient of the free energy. By com- 
paring with numerical derivatives, examples demonstrate that a reasonable 
accuracy is achieved. 

Keywords: analytical gradient, metals, free energy 
1. Introduction 

Today, analytical gradients are widely available in electronic structure 
codes. In the case of molecules, gradients with respect to the nuclear position 
are required, and in solids, in addition, gradients with respect to the cell 
parameters. Periodic systems often employ plane waves as basis functions, 
but local basis sets are also popular [J 0]. Local basis sets, usually atom 
centered, require the calculation of derivatives of the basis functions with 
respect to the nuclear positions, the Pulay forces [1, [ 3, III . This holds for the 
case of molecular and periodic 0, Q, 0, 0, M, Qllltll Q EE H E3, [l8, 



19l . |20|, |21| systems. Periodic systems have the feature that metallic ground 
states are a possible solution. Metallic systems are more difficult to treat than 
insulators, because the position of the Fermi energy has to be determined, and 
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the integration is only over a part of the Brillouin zone and thus more difficult 
than in the case of insulators. In the case of Hartree-Fock theory, there are 



further problems due to the vanishing density of states at the Fermi level [22 



and the slow decay of the density matrix at zero temperature (this is however 



less problematic at finite temperature where the decay is exponential 23|). 
This has motivated the use of a screened Coulomb operator for the exchange 
interaction 1241. For an overview of calculations for metals with Gaussian 



basis sets, see 25]. Some time ago, it had been argued that the gradient 



requires an extra term due to the shape of the Fermi surface [26|. This will 
be discussed in the present work, and it appears that this term is spurious. 
Numerical tests indicate that a reasonable accuracy can be achieved, and the 
analytical derivatives agree well with numerical derivatives of the free energy. 

2. Formalism 

2.1. Zero temperature 

The analytical gradients for periodic systems, on the Hartree-Fock level, 
were introduced by 0, 0] • A little later, an article suggested that an extra 
term should appear in the case of metals [26|, which will be reconsidered 
in the following. A notation similar to 0,0, 26 j is used, for the sake of 
simplicity. This corresponds to the case of one dimensional periodicity, but 
the argument can analogously be transferred to two and three dimensions. 



The notation is similar to the molecular case [27[, apart from the summation 
over the lattice vectors. 

The crystalline orbitals \l/n(r, k), with the band index n and the fc-point 
k are expanded in linear combinations of Bloch functions: 

# n (f, k) = J2 c M n(A;)^(r, k) (1) 

with 

k) = -J= exp(%'a)x^(r) (2) 
j 

where N is the number of unit cells in the macro-lattice, or equivalently 
the number of reducible fc-points, and xji( r ~) being a basis function (e.g. a 
Gaussian) in cell j. The overlap matrix element between orbital /i in cell 
and v in cell j is obtained as 



2 



S% = I xT(r)xl(rld 3 r 



(3) 



and its Fourier transform as 



S,u(k) = e*p(&ja)S% and S% = - ^ S^k) exp(-ikja) (4) 

3 k 

with the cell parameter a. Because of the orthonormality of the crystalline 
orbitals, it holds: 

c *nm(k)S 'fj, u (k)c un (k) = 5 mn (5) 
The total energy per primitive unit cell is expressed as in 0,0 , |26j as 

E = I E W + KiX + E ( NR ) ( 6 ) 

3,1*," 

with H®1 being the one-electron part of the Fock matrix element, F®i the 
corresponding Fock matrix element: 



pOj = TiOj , \^ plh ( 0j \ \hl 
fiv fiu ~ / j A 



At I \ t\J 



(7) 



with (^,11^) = (rlllut) being the two-electron integral: 



(Oj \ \hl\ 
\fiv II tA ) 



\ I xt{ri)xWi)r^^ 
2 J \r\-r 2 



Xr*(r2)x l x(r2)d 3 r l( i 3 r 2 
Xr*(r 2 )xi(r2)d 3 ri d s r 2 (8) 



n - r 2 



PI® is the corresponding density matrix element, and the nuclear repulsion 



energy is labelled as E(NR). Strictly speaking, some of the terms such as 
E(NR) are divergent for a periodic system, and a formulation based on e.g. 
the Ewald and related methods would be more suitable |28j,|29(. However, the 
main issue of the present paper can easiest be demonstrated with a notation 
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consistent with references 0, Q, |26| , and convergence issues of the Coulomb 
sums shall be ignored. The Hartree-Fock equations for periodic systems 
30, |3l| are: 



(9) 



with e n (k) being the eigenvalues. 
For metallic systems, the density matrix is expressed as in 26 



P 5 = I E ex P (ikja)c; n (k)c un (k)9(E F - e n (k)) (10) 

k,n 

k 

with the Fermi energy Ep and the Heaviside function 9. The factor 2 is 
due to the summation over the 2 spin states. Due to translational invariance, 
relations such as P^}„ = P„~ l hold. The derivative of the total energy with 
respect to a geometrical parameter |^ is then obtained as in p, |7|: 



dE_ = d^ffi pi0 , I V V P Ifc pjo dWri) n 1 n 

j,IJt,V j,fJ,,U h,l,T,\ 

The expression 



x exp (ikja)c; n (k)c„ n (k)9(E F - e n (k))e n (k) (12) 

fc,n 

corresponds to the energy weighted density matrix. 

In the following, the derivative of the 9 function shall be considered in 
more detail. When computing the gradient with respect to a geometrical 
parameter X, then the derivative term G due to the Heaviside function is 
obtained as 
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G = E 



dE dPff 89 (Ep - e n (k)) 



n ^dPffde(E F -e n (k)) 



E \ (TO 1 + F W + E p < 

m,a,p \ h,l,r,\ 



<9X 



A A a/9 I ItAJ 



J) pmO 



(13) 

06 (Ep - e n (k)) 



d6(E F - e n (k)) 



dX 



a/3 



pmO 

Ur l3a 



dO(E F - e n (k)) 



m,a,/3 



d6(E F - e n (k)) 



E iv F ^ E E ew(ikma)c* an (k)c Pn (k)5(E F - e n (k)) 



k n 



dEp de n (k)' 



dX 



dX 



^J2J2F a p(k)cUk)cp n (k)6(Ep - e n (k)) 



N 



k,n a,, 



E E c: n (k)S a p(k)cp n (k)6 n (k)5(Ep - e n (k)) 



dE F de n {k) 
dX dX 

8E F de n (k) 



k,n a,/3 



dX 



dX 



-Y,e n (k)6(E F -e n (k)) 



N 



k,n 



J2E F 5(Ep-e n (k)) 



k:n 



dE F de n (k) 
~dX~ ~ dX 

dEp de n (k) 



dX 



dX 



2 ^ d9{E F - e n (k)) 



N 



k,n 



dX 



Note that in reference (H % + F^L) appears instead of 2F^, and this 



appears to be incorrect (see also the related calculation in |27|). With the 



number of electrons in the unit cell no, it follows as in [26 



™o = E^ p 5 ( 14 ) 

2 . 

= E S °^m E exp (ikja)c^ n (k)c un (k)9(E F - e n (k)) 

k,n 

= Jj E E S^(kK n (k)c un (k)6(E F - e n (k)) = Q ( E * ~ e ^ 

[i,v k,n k,n 

and, as the particle number is fixed, ^ = 0, and therefore from equation 
[TBI G = is obtained: there is thus no extra term due to the step function, 
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and the same expression as for the case of insulators 0, 0] can be used for 
the derivatives with respect to geometrical parameters. 

2.2. Finite temperature 

An additional problem in the case of metals is the numerical integration of 
integrals over the occupied part of the Brillouin zone. This problem requires 
/c-point meshes as large as possible. A more efficient way is to apply a fi- 
nite temperature scheme. The calculation can then be theoretically based on 



finite temperature density functional theory 32]. The occupation numbers 



can be chosen e.g. accor ding to the Fermi function. Gaussian broadening is 



another popular scheme 33|, [3J, |35j . Further schemes (Lorentzian broaden- 



ing, a step function) had been discussed in [36| . The Fermi function has the 
advantage that the computed free energy has a direct physical meaning, as 
it contains the electronic contribution to the free energy; contributions due 
to e.g. phonons are however missing (see, e.g. ji?}). The Fermi function is 
defined as 



fk,r 



1 



(15) 



l + exp((e n (£;) -E F )/k B T) 

with the Boltzmann constant k%- A small finite temperature can be 
introduced, so that the density matrix becomes 



pOj 



N 



exp {-ikja)c* vn {k)c^ n {k)f Kr , 



(16) 



k,n 



and 



Pfiu(k) — 2 y ] C un (k)Cfj,n(k) fk,n 



(17) 



Compared to equation [TO], the Heaviside function was replaced with the 
Fermi function. At zero temperature, the equations agree. The zero temper- 
ature energy can subsequently be approximated by [38| 



E(0) = h(E(T)+F(T)) 



(18) 



with the entropy 
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S(T) 



2k f 



N 



J^(/*,n ^ fk,n + (1 - fk,n) - /fc, n )) 



(19) 



k.n 



and the free energy 



F(T) = E(T) - TS'(T) (20) 

F(T) and -E'(T) are similar at low temperature, and the error should be 
relatively small when using F(T) instead of E(T). As was pointed out later 



40, 41 



analytical gradients are, for the case of an occupancy according 
to the Fermi function, consistent with the free energy F(T). This can be 
seen by computing the additional terms due to the entropy: 



dS(T) _ 2k B T v df k , n f k ,n (21) 

OX ' N j£ dX n i-f k , n 1 ' 

2 V—. df k n 2 ^ df k n 

k.n k,n 

Here, it was exploited that 4 n f k) n = hi analogy to equation [141 

and thus the derivative jfYlkn'Bx^^F = ®- Another term is due to the 
derivative of the density matrix. 

This leads now to an additional term: 



E dE \^ dPf a df k>n 
QpjO Z^Qf qx { > 

j,a,P 13a k,n Jlt ' n 

j,a,p P a k,n ' 

= Jf E eXP ( i ^ a ) F S C an( fc )C/3n(fc)^ !i 

j,a,/3 k,n 

= ^^ c ln{ k )S a p{k)c Pn {k)t n {k)^^- = ^X^^l^ 

j,a,P k,n 

But this term is just equivalent to the entropy term in equation [211 with 
opposite sign. As a whole, for the derivatives of the free energy with respect 
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to a geometrical parameter X, the two terms containing derivatives of the 

9fk,n 

ax 



occupation number cancel, and the expression is: 



dF dE-TS dH% nm , 1 ^ ^ ^ ^ d(%\\¥xJ 



dX dX ^ dX 2 ^ ^ AT dX 

h,l,r,\ 

- E If E I ex P (rtja)c; n (k)c un (k)f k , n e n (k) + (23) 



This can be viewed as a generalization of the result in section I2.1[ with 
the 9 function being replaced with the Fermi function. At zero temperature, 
this reduces to the 9 function, and the entropy becomes zero. These argu- 
ments hold similarly for the case of higher dimensions or the case of density 
functional theory. 

For higher temperatures T, the forces and the derivative of the total 
energy deviate stronger, and a suggestion was made to remedy this, in order 



to obtain the derivative of the total energy, and not of the free energy [42 



3. Examples 

In the following, some examples demonstrate the accuracy of the gra- 
dients. The calculations were done with the present CRYSTAL09 release 
43, The examples aim at documenting the accuracy of the gradient, by 



comparing the analytical and numerical gradient, at the level of Hartree- 
Fock and density functional theory, for the gradient with respect to the cell 
parameter, and with respect to the nuclear position. 

First, for Cu bulk, the analytical and numerical gradient with respect to 
the cell parameter are compared in table [TJ This is done on the Hartree-Fock 



and density functional level. The basis sets from reference |44| were used. A 
fc-point mesh with 16 x 16 x 16 points was used. Smearing temperatures in 
the range from 0.001 to 0.05 E^ were chosen. Technically, in the input, 
a hybrid functional consisting of nothing but 100% Fock exchange was de- 
fined, in order to perform the Hartree-Fock calculation at finite temperature. 
When comparing numerical and analytical derivatives, then the obtained ac- 
curacy for the derivative of the free energy — ^ is similar to the one for 



insulators, see 15, 16, 17, 18, 191. Note that in addition, the numerical noise 



is in general larger in the case of metals, and therefore, also the energies and 
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Table 1: The derivative of the total energy and the free energy, in hartree/bohr (Eh/ao), 
with respect to the cell parameter a, analytical and numerical, on the Hartree-Fock and 
density functional (LDA) level. 



smearing temperature — |p (numerical) — |p (numerical) — |p (analytical) 
(E h ) a & a & a & 

Hartree-Fock (at a = 5 A) 



0.001 -0.0316 -0.0316 -0.0314 

0.01 -0.0317 -0.0315 -0.0313 

0.03 -0.0328 -0.0305 -0.0303 

0.05 -0.0352 -0.0276 -0.0280 

LDA (at a = 3.4 A) 

0.001 0.0315 0.0315 0.0317 

0.01 0.0310 0.0319 0.0320 

0.03 0.0212 0.0390 0.0393 

0.05 0.0098 0.0540 0.0542 



their numerical derivatives carry larger noise. The agreement between ana- 
lytical and numerical derivative of the free energy is similar for all smearing 
temperatures. 

The derivative of the energy with respect to the cell parameter agrees 
reasonably well at low temperatures, but deviates strongly at high smearing 
temperatures, as expected, as the energy and the free energy deviate more 
and more at higher temperature. The free energy and its derivative with 
respect to the cell parameter are also visualized in figure HJ where a smear- 
ing temperature of 0.001 Eh was employed. Again, the agreement between 
numerical and analytical derivative is very good. 

As an example for the gradient with respect to nuclear positions, the ad- 
sorbate system Cu(lll)(v3x \/3)R30 o -Cl is considered, with chlorine sitting 
on the hep (hexagonal close packed) site. The basis sets are as in 44j, and 



16 x 16 fc-points together with a smearing temperature of 0.001 Eh is used. 
The free energy and its derivative with respect to the z-component of the 
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chlorine atom are computed analytically and numerically. The results are vi- 
sualized in figure HJ and the numerical and analytical derivatives agree well. 
The computed equilibrium position corresponds to a hight of 1.85 A above 
the topmost Cu layer, in reasonable agreement with the earlier calculation 



441 ]: in the earlier calculation, a generalized gradient functional had been 
employed and a hight of 1.90 A had been obtained. The present calculation 
gives a slightly shorter bond length which is a usual feature of the local den- 
sity approximation (LDA), as compared to gradient corrected functionals. 
Note that no gradients had been used in the earlier work [44J, and the geom- 



etry had been determined by iteratively optimizing the various geometrical 
parameters, by employing the total energy only. 



4. Conclusion 

Derivatives of the total and free energy of periodic systems with respect 
to geometrical parameters were studied theoretically, in the case of metallic 
systems. In the case of metals, numerical integration is often facilitated by 
introducing an artificial temperature and by an occupancy according to e.g. 
the Fermi function. At zero temperature, the theory of the derivatives does 
not require an additional term compared to the case of insulators. At finite 
temperature, when the occupancy is according to the Fermi function, then a 
similar expression for the derivative can be employed, which is however only 
consistent with the free energy. Therefore, numerical derivatives of the free 
energy agree reasonably well with analytical derivatives, and consequently, 
numerical derivatives of the total energy deviate more and more with in- 
creasing temperature. This holds for the case of Hartree-Fock or density 
functional theory. Numerical examples demonstrate the accuracy which is 
achieved with the implementation in the CRYSTAL code. 
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Figure 1: Left: Free energy for Cu bulk; crosses refer to computed points, the full line is 
a fit through the points. Right: Analytical (crosses) and numerical derivative (full line) 
with respect to the cell parameter for Cu bulk. The numerical derivative is obtained as a 
derivative of the fit of the energy expression in the left figure. A smearing temperature of 
0.001 Eh was applied. 



Cu: free energy Cu: derivative of the free energy with respect to the cell parameter 

LDA LDA 




lattice constant (A) lattice constant (A) 
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Figure 2: Left: Free energy for Cl/Cu(lll); crosses refer to computed points, the full line 
is a fit through the points. Right: Analytical (crosses) and numerical derivative (full line) 
with respect to the z-position of the CI atom. The numerical derivative is obtained as a 
derivative of the fit of the energy expression in the left figure. A smearing temperature of 
0.001 Eh was applied. 




Cl/Cu(l 11): force on CI atom, z component 
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